#ifndef USE_mpi
void checkOMP(){

#ifdef USE_omp

    int nthreads = 0; 
    int th_id;
#pragma omp parallel private(th_id) shared(nthreads)
    {
	th_id = omp_get_thread_num();
#pragma omp critical
	{
	    std::cout << "Hello World from thread " << th_id << '\n';
	}
#pragma omp barrier

#pragma omp master
	{
	    nthreads = omp_get_num_threads();

	    std::cout << "There are " << nthreads << " threads" << '\n';
	}
    }

#endif 
  
}
#endif 

void printvec(double x[], int n){
	for (int i=0; i<n; i++){
		printf("\n i %d, %f ", i, x[i] );	
	}	
}
void printvec(gsl_vector* x, int n){
	for (int i=0; i<n; i++){
		printf("\n i %d, %f ", i, gsl_vector_get(x, i) );	
	}	
}

void printvec(const char *fname, double x[], int n){

	FILE * outfile; outfile  = fopen (fname, "w"); 
	for (int i=0; i<n; i++){
		fprintf(outfile, "%.8f \n", x[i] );	
	}	
	fclose(outfile); 	
}

void printvec(const char *fname, int x[], int n){

	FILE * outfile; outfile  = fopen (fname, "w"); 
	for (int i=0; i<n; i++){
		fprintf(outfile, "%d \n", x[i] );	
	}	
	fclose(outfile); 	
}

void printvec(const char *fname, gsl_vector* x, int n){

	FILE * outfile; outfile  = fopen (fname, "w"); 
	for (int i=0; i<n; i++){
		fprintf(outfile, "%.8f \n", gsl_vector_get(x, i) );	
	}	
	fclose(outfile); 	
}
void printvec(const char *fname, const gsl_vector* x, int n){

	FILE * outfile; outfile  = fopen (fname, "w"); 
	for (int i=0; i<n; i++){
		fprintf(outfile, "%.8f \n", gsl_vector_get(x, i) );	
	}	
	fclose(outfile); 	
}


// Dynamic array allocation
// prevent orphan pointer problem;
int **my_allocate_int(int rows, int cols){
	int **array_my_allocate_int; 
    array_my_allocate_int = new int* [rows]; 
    for(int i = 0; i < rows; i++) array_my_allocate_int[i] = new int [cols];
    return array_my_allocate_int;
}
double **my_allocate_db(int rows, int cols){
	double **array_my_allocate_db; 
	array_my_allocate_db = new double* [rows];
    for(int i = 0; i < rows; i++){
		array_my_allocate_db[i] = new double [cols];
		if (array_my_allocate_db[i] == NULL) std::cout << "Could not allocate memory" << std::endl;
	}
    return array_my_allocate_db;
}
void my_free(double** ptr, int rows){
    for(int i = 0; i < rows; i++)
		delete [] ptr[i];
}
void my_free(int** ptr, int rows){
    for(int i = 0; i < rows; i++)
		delete [] ptr[i];
}


int my_sum(const int N, int x[]){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += x[i];
    return my_sum; 
}
double my_sum(const int N, double x[]){
    double my_sum = 0.0;  
    for (int i=0; i< N; i++)	my_sum += x[i];
    return my_sum; 
}

//--- my_count() returns the number of x that is within [xmin, xmax] 
int my_count(const int N, int x[], int xmin, int xmax){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] >= xmin ) * ( x[i] <= xmax);
    return my_sum; 
}
int my_count(const int N, double x[], double xmin, double xmax){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] >= xmin ) * ( x[i] <= xmax);
    return my_sum; 
}


int my_countMore(const int N, double x[], double xmin){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] > xmin );
    return my_sum; 
}
int my_countMore(const int N, int x[], int xmin){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] > xmin );
    return my_sum; 
}
int my_countLess(const int N, double x[], double xmax){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] < xmax);
    return my_sum; 
}
int my_countLess(const int N, int x[], int xmax){
    int my_sum = 0;  
    for (int i=0; i< N; i++)	my_sum += ( x[i] < xmax);
    return my_sum; 
}



double my_pct(const int N, const double x0[], const int p){
//... it is crucially important to keep x0[] as constant 
//... sort() will change the order of input x[]
   	
   	if (p >100 or p<0){ printf("\n Errors in my_pct: p %2d \n\n existing \n\n", p);  exit(1); }
   	
   	double x[N]; 
   	for (int i=0; i< N; i++) x[i] = x0[i]; 
	   	
	sort (x, x + N); 
	   	
   	double pct = 0.0; 
   	
   	int np = N*p/100; 
	    
    if ( N*p % 100 ==0 )	{
      pct = (x[ np - 1] + x[ np]) / 2.0;
  	} else 	{
      pct = x[ np];
  	}
	
//	cout << "Sorted Array looks like this." << endl;
//  for (int i = 0; i != N; ++i) cout << x[i] << " ";
//	cout << "\n\t median = " << median << endl;
	        
    return pct; 
}


double my_median(const int N, const double x0[]){
//... it is crucially important to keep x0[] as constant 
//... sort() will change the order of input x[]
   
   	double median = 0.0; 
   	
   	double x[N]; 
   	for (int i=0; i< N; i++) x[i] = x0[i]; 
	   	
	sort (x, x + N);    	
	
    if (N % 2 == 0)	{
      median = (x[ N / 2 - 1] + x[ N / 2]) / 2.0;
  	} else 	{
      median = x[ N/ 2];
  	}
	
//	cout << "Sorted Array looks like this." << endl;
//  for (int i = 0; i != N; ++i) cout << x[i] << " ";
//	cout << "\n\t median = " << median << endl;
	        
    return median; 
}


double my_mean(const int N, double x[]){
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += x[i]; 
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(N+0.0); 
}
double my_mean(const int N, int x[]){
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += x[i]*1.0;
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(N+0.0); 
}
double my_var(const int N, double x[]){

    double meanx = my_mean(N, x);
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += (x[i]-meanx)*(x[i]-meanx);
    return ( my_sum/(N-1) ); 
}
double my_var(const int N, int x[]){
	
    double meanx = my_mean(N, x);
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += (x[i]-meanx)*(x[i]-meanx);
    return ( my_sum/(N-1) ); 
}
double my_sd(const int N, double x[]){

    double meanx = my_mean(N, x);
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += (x[i]-meanx)*(x[i]-meanx);
    return sqrt( my_sum/(N-1) ); 
}
double my_sd(const int N, int x[]){
	
    double meanx = my_mean(N, x);
    double my_sum = 0.0; 
    for (int i=0; i< N; i++)	my_sum += (x[i]-meanx)*(x[i]-meanx);
    return sqrt( my_sum/(N-1) ); 
}

double my_dotproduct(const int N, double x1[], double x2[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += x1[i]*x2[i]; 
	return s; 
}
double my_dotproduct(const int N, int x1[], double x2[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += 1.0*x1[i]*x2[i]; 
	return s; 
}
double my_dotproduct(const int N, double x1[], int x2[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += 1.0*x1[i]*x2[i]; 
	return s; 
}
double my_dotproduct(const int N, int x1[], int x2[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += 1.0*x1[i]*x2[i]; 
	return s; 
}
double my_dotproduct3(const int N, int x1[], int x2[], int x3[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += 1.0*x1[i]*x2[i]*x3[i]; 
	return s; 
}
double my_dotproduct3(const int N, double x1[], int x2[], int x3[]){
	double s=0.0;
    for (int i=0; i< N; i++) s += 1.0*x1[i]*x2[i]*x3[i]; 
	return s; 
}

double my_skewness(const int N, double avg, double sd, double v[]){
	
	double skewness = 0.0;
    for (int i = 0; i < N; i++){ skewness += (v[i] - avg)*(v[i] - avg)*(v[i] - avg); }
    skewness = skewness/(N * sd * sd * sd);
	return skewness; 	
}

double my_mean_conditional(const int N, double x[], int cat[], int z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] == z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_mean_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditional(const int N, int x[], int cat[], int z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] == z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_mean_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditional(const int N, double x[], int cat[], int zlow, int zhigh){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] >= zlow && cat[i] <= zhigh){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_mean_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditional(const int N, int x[], int cat[], int zlow, int zhigh){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] >= zlow && cat[i] <= zhigh){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_mean_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}

double my_meanlog_conditional(const int N, double x[], int cat[], int zlow, int zhigh){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] >= zlow && cat[i] <= zhigh && x[i] > 1.0e-30){
    	my_sum += log(x[i]); 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_meanlog_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_meanlog_conditional(const int N, int x[], int cat[], int zlow, int zhigh){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] >= zlow && cat[i] <= zhigh && x[i] > 1.0e-30){
    	my_sum += log(x[i]); 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n Errors in my_meanlog_conditional: sumz <= 1 \n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}

double my_mean_conditionalLess(const int N, double x[], double cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] < z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalLess(const int N, int x[], double cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] < z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalLess(const int N, double x[], int cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] < z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalLess(const int N, int x[], int cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] < z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}




double my_mean_conditionalMore(const int N, double x[], double cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] > z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalMore(const int N, int x[], double cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] > z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalMore(const int N, double x[], int cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] > z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}
double my_mean_conditionalMore(const int N, int x[], int cat[], double z0){
    double my_sum = 0.0; double my_sumz = 0.0; 
    for (int i=0; i< N; i++){
    	if (cat[i] > z0){
    	my_sum += x[i]; 
    	my_sumz ++; 
    	}
	}
	if (my_sumz <= 1){ printf("\n\n Errors in my_mean_conditionalLess: sumz <= 1 \n\n exiting ... \n\n"); }
//    printf("\n my_sum %f, mean %f \n", my_sum, my_sum/(N+0.0)); 
    return my_sum/(my_sumz+0.0); 
}



double my_gini(const int N, const double x[]){
// ... calculating Gini coefficient, might be above 1 if there are negative values. 
	
	double gini = 0.0; double sumx = 0.0; 
	
	for (int i1=0; i1<N; i1++){
		for (int i2=0; i2<N; i2++){	
			gini += fabs(x[i1] - x[i2]); 
		}
		sumx += x[i1]; 
	}
		
	gini /= (2.0*sumx*N); 
	
	return gini; 
}

double my_gini_positive(const int N, const double x[]){
// ... calculating Gini coefficient for all positive values
    
    const double smallpos = 1.0e-8;
		
	double gini = 0.0; double sumx = 0.0; 
	
	for (int i1=0; i1<N; i1++){
		for (int i2=0; i2<N; i2++){	
		
			if (x[i1] > smallpos && x[i2] > smallpos){ 
			
			gini += fabs(x[i1] - x[i2]); 
			
			}
		}
		sumx += x[i1]; 
	}
		
	gini /= (2.0*sumx*N); 
	
	return gini; 	

}

double my_varlog(const int N, int x[]){
//  very small numbers are ignored 
    
    const double smallpos = 0.0001;
    
    int npos = 0; 
    for (int i=0; i< N; i++){
        if (x[i] > smallpos) npos++;
    }
    
    double logx[npos];
    
    int inpos = 0; 
    for (int i=0; i< N; i++){
    
        if (x[i] > smallpos){ 
        	logx[inpos] = log(x[i]);
        	inpos++; 
        }
    }
    
    double meanlogx = my_mean(npos, logx);
    
    double my_sum = 0.0;
    for (int i=0; i< npos; i++)	my_sum += (logx[i]-meanlogx)*(logx[i]-meanlogx);
    
    return ( my_sum/(npos-1) );
}


double my_varlog(const int N, double x[]){
//  very small numbers are ignored 
    
    const double smallpos = 0.0001;
    
    int npos = 0; 
    for (int i=0; i< N; i++){
        if (x[i] > smallpos) npos++;
    }
        
    double logx[npos];
    
    int inpos = 0; 
    for (int i=0; i< N; i++){
    
        if (x[i] > smallpos){ 
        	logx[inpos] = log(x[i]);
        	inpos++; 
        }
    }
    
    double meanlogx = my_mean(npos, logx);
    
    double my_sum = 0.0;
    for (int i=0; i< npos; i++)	my_sum += (logx[i]-meanlogx)*(logx[i]-meanlogx);
    
    return ( my_sum/(npos-1) );
}


double my_varlog1(const int N, double x[]){
//  very small numbers are ignored 
        
    double logx[N];
    
    for (int i=0; i< N; i++){
       	logx[i] = log(x[i]+1.0);
    }
    
    double meanlogx = my_mean(N, logx);
    
    double my_sum = 0.0;
    for (int i=0; i< N; i++)	my_sum += (logx[i]-meanlogx)*(logx[i]-meanlogx);
    
    return ( my_sum/(N-1) );
}

double *gen_Grid(const int n, double c_min, double c_max, double cc){
//  Generate a grid with length n, between [c_min, c_max] with curvature cc

	#ifdef DEBUG
		if (n<2) { printf("\n\n Errors in gen_Grid(): n %d < 2 \n\n exiting \n", n); exit(1); }
	#endif
		
	double *grid; grid = new double [n]; 	
		
	for (int i = 0 ; i < n; i++){
		
		double frac = pow( i /(n-1.0), cc ); 
		
		grid[i] = (c_max - c_min) * frac + c_min;  
	
//		printf("\n i %d, grid %f, %f ", i, grid[i], frac); 
			
	}
	
	return grid; 
}


double *gen_Grid2(const int n_left, const int n_right, double c_min, double c_max, double c_central, double cc){
//  Generate a grid with length n, between [c_min, c_max], concentrated around c_central with curvature cc
	
	if (c_central < c_min || c_central > c_max || cc < 1.0){
		printf("\n Errors in gen_Grid(): n_left %d, n_right %d, c_min %f, c_max %f, c_central %f, cc %f \n ... exiting ...\n", 
		n_left, n_right, c_min, c_max, c_central, cc);
		exit(1);  
	}
	
	const int n = n_left + n_right; 
	double *grid; grid = new double [n]; 	
	
	for (int i = 0 ; i < n_left ; i++){
		
		double frac = pow( (i+1.0) /n_left, cc ); 
		
		int igrid = n_left-1- i;  
		
		grid[igrid] = (c_min - c_central) * frac  + c_central;  
	
//		printf("\n i %d, grid[%d] %f, %f ", i, igrid, grid[igrid], frac); 
	
	}
	
	
	for (int i = 0 ; i < n_right ; i++){
		
		double frac = pow( i /(n_right-1.0), cc ); 
		
		int igrid = i + n_left ; 
		
		grid[igrid] = (c_max - c_central) * frac + c_central;  
	
//		printf("\n i %d, grid[%d] %f, %f ", i, igrid, grid[igrid], frac); 
			
	}
	
	return grid; 
}


// used in GS filter of standard error calculation. 
int my_index_stderror(const int z, const int np){
//  seek shift index for given case nl, nr, m (see savgol).
	int ic = 0; 
	if (z<0) ic = -z; 
	if (z>0) ic = np - z; 
	return ic; 
}

//!!!!! The following two functions: olsReg() and olsRsquare(), depend on GSL library!!!!!

double olsReg(const int n, const int nx, const gsl_vector *y, const gsl_matrix *X, gsl_vector *coef){
// This function calculate the coefficients of a multiple linear regression model using least squares.
	
    gsl_matrix *ols_cov  	= gsl_matrix_alloc(nx, nx); 

    double chisq;	

    // Allocate a workspace for polynomial fitting: 
    gsl_multifit_linear_workspace * work =  gsl_multifit_linear_alloc (n, nx);	
	
    // Initialization
    gsl_vector_set_zero(coef); 

    // Get coefficients
    gsl_multifit_linear(X, y, coef, ols_cov,  &chisq, work); 
	
	
#ifdef DEBUG
    // Check matrix inverse of the OLS X matrix, 
    // Singular value decomposition (SVD), A = U S V', the presence of a zero in S indicates that the matrix is singular 
    // check the determinant of (X'X), if the determinant of X'X is close to 0, then (X'X) is near singular
	
    // X = U S V^T
    gsl_matrix * U    = gsl_matrix_alloc(n, nx); 
    gsl_matrix * V    = gsl_matrix_alloc(nx, nx); 
    gsl_vector * S    = gsl_vector_alloc(nx); 
    gsl_vector * W = gsl_vector_alloc(nx); 
	
    // U = X; 
    gsl_matrix_memcpy( U, X) ; 
    gsl_linalg_SV_decomp (U,  V,  S, W); 
	
    // only print the smallest root
    double temp = gsl_vector_get(S, nx-1); 
    if (temp < 1.0 ){
    printf("\n \t Warning: check singularity of ols regression matrix: S[i] = %f ", gsl_vector_get(S, nx-1) ); 
	}
	
    gsl_matrix_free(U); 
    gsl_matrix_free(V);
    gsl_vector_free(S); 
    gsl_vector_free(W);
	
#endif
	
    // free memory from the multi-parameter fitting workspace
    gsl_multifit_linear_free (work); gsl_matrix_free(ols_cov);
    
    return chisq; 
}


double olsRsquare(const int n, const gsl_vector *y, double chisq){

    double ymean, tss, rsquare; 
		
    ymean = 0.0; 
    for (int i=0; i<n; i++) {
	ymean += gsl_vector_get(y,i); 
    }
    ymean /= (1.0*n); 
		
		
    tss = 0.0; 	
    for (int i=0; i<n; i++){
		
	tss += (gsl_vector_get(y,i)-ymean)*(gsl_vector_get(y,i)-ymean); // TSS: total sum of squares
		
    }

    rsquare = 1.0-chisq/tss; 
	
    return rsquare; 
}

double my_min(const int N, const double x[]){
   
   	double y = x[0]; 
	   	
   	for (int i=0; i< N; i++){
   		if ( y > x[i] ) y = x[i]; 
   	}
   	        
    return y; 
}
int my_min(const int N, const int x[]){
   
   	int y = x[0]; 
	   	
   	for (int i=0; i< N; i++){
   		if ( y > x[i] ) y = x[i]; 
   	}
   	        
    return y; 
}
double my_max(const int N, const double x[]){
   
   	double y = x[0]; 
	   	
   	for (int i=0; i< N; i++){
   		if ( y < x[i] ) y = x[i]; 
   	}
   	        
    return y; 
}
int my_max(const int N, const int x[]){
   
   	int y = x[0]; 
	   	
   	for (int i=0; i< N; i++){
   		if ( y < x[i] ) y = x[i]; 
   	}
   	        
    return y; 
}
double my_wmean_omp(const int N, int x[], double wgt[]){
    
    double my_sum = 0.0; 
    double twgt = 0.0; 
        
	#ifdef USE_omp
	#pragma omp barrier
	#pragma omp flush
				
	#pragma omp parallel for reduction(+: my_sum, twgt)
	
	#endif
	
    for (int i=0; i< N; i++){
    	my_sum += x[i] * wgt[i];
    	twgt   += wgt[i];
    }

	#ifdef USE_omp
	#pragma omp barrier
	#pragma omp flush
	#endif
	    
    return ( my_sum/twgt ); 
}
double my_wmean_omp(const int N, double x[], double wgt[]){
    
    double my_sum = 0.0; 
    double twgt = 0.0; 
        
	#ifdef USE_omp
	#pragma omp barrier
	#pragma omp flush
				
	#pragma omp parallel for reduction(+: my_sum, twgt)
	
	#endif
	
    for (int i=0; i< N; i++){
    	my_sum += x[i] * wgt[i];
    	twgt   += wgt[i];
    }

	#ifdef USE_omp
	#pragma omp barrier
	#pragma omp flush
	#endif
	    
    return ( my_sum/twgt ); 
}

double my_wmean_omp_test(){
	
	int n = 100; 
	double x[n], wgt[n]; 
	
	for (int i=0; i<n; i++){
		
		x[i] = i/n + i; 
		wgt[i] = i*i/n; 
	}

	double wm = my_wmean_omp(n, x, wgt); 
	
	double wm1; 
	
    double my_sum = 0.0; 
    double twgt = 0.0; 
    for (int i=0; i< n; i++){
    	my_sum += x[i] * wgt[i];
    	twgt   += wgt[i];
    }
	wm1 = my_sum/twgt; 
	
	printf("\n\n test my_wmean_omp() \n \t wm %f, wm1 %f \n exiting ...\n", wm, wm1); exit(1);  


}

double cond_mean_cutoff(const double cutoff,const double sigma_eta){
// ... E_eta: mean of the truncated normal, "Inverse Mills Ratio" 
	
	double e_eta = 0.0; 
	 	
 	double z = cutoff / sigma_eta; 
 	
 	
 	double pdf = gsl_ran_ugaussian_pdf ( z );  // standard normal density function
 											   // equal to sigma_eta * gsl_ran_gaussian_pdf ( cutoff, sigma_eta ); 
 											   
 	double cdf = gsl_cdf_ugaussian_P ( z );    // standard normal cumulative distribution function


// 	      double pi = 3.141592653589793;
//        E_eta=-dexp(-0.50*(cutoff/sigma_eta)**2.0)/sigma_eta/(2.0*ppi)**0.50
//        E_eta=E_eta/normcdf(dble(cutoff/sigma_eta), .FALSE.)*sigma_eta
	
	e_eta = - ( pdf / cdf ) * sigma_eta; 
	
	return e_eta; 
}
